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We present a method, based on characterizing efficiency fluctuations, to asses the performance 
of nanoscale thermoelectric junctions. This method accounts for effects typically arising in small 
junctions, namely, stochasticity in the junction’s performance, quantum effects, and nonequilibrium 
features preventing a linear response analysis. It is based on a nonequilibrium Green’s function 
(NEGF) approach, which we use to derive the full counting statistics (FCS) for heat and work, and 
which in turn allows us to calculate the statistical properties of efficiency fluctuations. We simulate 
the latter for a variety of simple models where our method is exact. By analyzing the discrepancies 
with the semi-classical prediction of a quantum master equation (QME) approach, we emphasize the 
quantum nature of efficiency fluctuations for realistic junction parameters. We finally propose an 
approximate Gaussian method to express efficiency fluctuations in terms of nonequilibrium currents 
and noises which are experimentally measurable in molecular junctions. 


PACS numbers: 05.70.Ln 85.65.+h 85.80.Fi 84.60.Rb 

I. INTRODUCTION 

The development of thermoelectric materials is at the 
forefront of the research related to energy conversion and 
storage. While research on thermoelectricity in bulk ma¬ 
terials goes back to the middle of the last century jl|, 
measurements at the nanoscale (and in particular, stud¬ 
ies of thermoelectricity in molecular junctions) were only 
reported recently Q. The small size of the junctions 
gives rise to new physical phenomena, not accessible at 
the macroscopic level, and which are considered promis¬ 
ing for reaching more effective energy conversion. The 
thermoelectric properties of nanoscale junctions have in¬ 
deed received a lot of attention in the last years, both 
experimentally Ml and theoretically ITbl [2!)t . 

Experimental studies on thermoelectricity in nanoscale 
junctions make use of the macroscopic theory of thermo¬ 
electricity to asses the junction’s performance. The latter 
is characterized by the figure of merit, a quantity exclu¬ 
sively defined in terms of linear response transport coef¬ 
ficients and thus ill-defined out of nonequilibrium. While 
the linear theory is reasonable in bulk material, it fails 
in small thermoelectric junctions which can operate in 
the nonlinear regime (for instance in the resonant tun¬ 
neling regime). This fact motivated a number of studies 
to consider the macroscopic efficiency of the junction as 
an alternative to the figure of merit to characterizes the 
performance of the junction [do! Pidl . The macroscopic ef¬ 
ficiency is the traditional thermodynamic efficiency of a 
heat engine defined as the fraction of average power out¬ 
put extracted from the heat arising from the hot source. 

It is well defined far from equilibrium and upper bounded 
by the Carnot efficiency. 

The nonequilibrium features of the junction are not the 
only characteristic to be accounted for at the nanoscale. 


Due to the small size of the system, thermal fluctuations 
will play a much more import role than in bulk samples, 
resulting in a high variability in the junction’s perfor¬ 
mance. This variability requires a statistical characteri¬ 
zation of the energy conversion which can be performed 
using the methods of stochastic thermodynamics [44H47j . 
Such studies have been recently done for small classical 
energy converters 0-Ell • The main idea is to define the 
efficiency along a single realization of the operating de¬ 
vice and to develop techniques to study its fluctuations. 
Experimental studies of efficiency fluctuation have been 
very recently performed in Ref. |54| . 

The third central feature of small thermoelectric junc¬ 
tions which needs to be accounted for are quantum ef¬ 
fects. Indeed, quantum coherences can significantly af¬ 
fect charge and energy transfers in molecular junctions 
as discussed theoretically in Refs [55l - l60| and shown ex¬ 
perimentally in Refs. (6lM66j| . 

In this paper we provide a general method to study 
the performance of nanoscale thermoelectric junctions 
based on efficiency fluctuations. This methods accounts 
for the three key features characterizing small junctions, 
namely, the variability in performance due to fluctua¬ 
tions, operation modes arbitrary far from equilibrium, 
and quantum effects. It is based on the joint energy and 
particle full counting statistics (FCS) which we calculate 
within the nonequilibrium Green’s functions (NEGF) for¬ 
malism. While the quantum FCS of particle currents 
(e.g. electrons) in junctions is well developed (67l - [78| . 
that of energy was mostly limited to the quantum mas¬ 
ter equatio n (Q ME) approach i 7ll f79M8^| with its known 
limitations [761, [85[ [86|. By numerically calculating effi¬ 
ciency fluctuations for a set of simple models and com¬ 
paring our NEGF results with those obtained using a 
QME approach, we identify the regimes where efficiency 
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fluctuations display truly quantum features. Moreover, 
we propose an approximate Gaussian scheme enabling 
to estimate efficiency fluctuations solely based on exper¬ 
imentally measurable quantities in molecular junctions 
[87l - f92l ] . namely the nonequilibrium energy and matter 
currents and noises. 

The structure of the paper is the following. After intro¬ 
ducing the FCS of energy, work, and heat within NEGF 
in Section IE we consider efficiency fluctuations in Sec¬ 
tion IIIII Iii Section IIVI we numerically evaluate efficiency 
fluctuations for various models, compare our results with 
the QME approach and describe the approximate scheme 
to estimate efficiency fluctuations experimentally. We 
summarize our findings in section Ivl 


II. FCS OF PARTICLE AND ENERGY FLUXES 


The particle FCS for a single level strongly coupled to 
Fermi reservoirs was derived in Ref. [blj and generalized 
to a multilevel interacting system in Ref. [75]. Later, the 
methodology was applied to describe inelastic transport 
in junctions in Ref. [77|, where the role of quantum co¬ 
herence on the FCS was discussed. 

Here we extend the methodology to count particles 
and energy in a system strongly coupled to its reservoirs. 
Similar to the particle FCS [7l[, the treatment starts 
by dressing the evolution operator, with particle, 

7 ^, and energy, 7 ^, counting fields at interface K of the 
junction 

U 7 (t,f') = e ~ i ^ K e- i ' y ^ 6K U(t,t , )e +i ' y ^ K e +i ' y ^ flK 

(!) 

Note that [ Nk\Hk] = 0. The counting fields depend on 
the Keldysh contour branch (see Fig. [lji) 


7 k 


d-A|(-/2 at — e ( +Aj|/2 at 
-A£/2 at + 7k = |-A§/2 at + 


Here — and + are the time ordered and anti-time ordered 
branches of the contour, respectively. 

Following the procedure outlined in Refs. m, m , at 
steady state we get the following expression for deriva¬ 
tives of the cumulant generating function, S = — i(tf — ti ) 
U (here U is the adiabatic potential), in the counting 
fields (Af = P, E ) 

= — j -^OmIk(E) (3) 

where Om = 1 ( E ) for M = P ( E ), and 

Ik{E) =Tr{s <(E)e i ^ +EX ^ ) G>(E) (4) 

-G <(E)S>(K)e" i(A « + - EA ^ ) | 

is the energy resolved dressed particle current at inter¬ 
face K, Tr{...} is the trace over the system subspace, 




FIG. 1: (Color online) (a) Sketch of the counting field A dress¬ 
ing of the Keldysh contour forward (—) and backward (+) 
branches, (b) Sketch of the a nano-thermoelectric junction 
consisting of a molecule M embedded between two contacts 
L and R with Tl < Tr and fiL > Hr- 


and > is the lesser (greater) projections of the Green 
function obtained from a counting field dressed version 
of the Dyson equation (see e.g. Ref. for details). 

While the expression for the derivatives of the adia¬ 
batic potential in the counting fields can be easily formu¬ 
lated in terms of the field-dressed Green functions and 
self-energies (see Eq. (HI) above or Ref. (77j for a detailed 
discussion), the corresponding expression for the adia¬ 
batic potential itself is more complicated. An explicit 
expression for the adiabatic potential within the NEGF 
based particle FCS for a single noninteracting level was 
derived in Ref. [Hj]. Exact results for the particle FCS 
for one-dimensional tight-binding junction models were 
presented in Ref ff^}. Here we consider the case of a 
multilevel noninteracting system. In particular, we show 
that for a single level coupled to its reservoirs (and pos¬ 
sibly also to other levels with the latter not coupled to 
reservoirs) or for a multilevel system coupled to reservoirs 
through single molecular orbitals, the explicit expression 
for the adiabatic potential in the presence of both particle 
and energy counting fields is (see Appendix for details) 

^({A}) = * J gln(l + T{E) 

x {/l(E)[1 - f R (E)][e + < x ^- x ^ +E ^- x ^ - 1] (5) 

+f R (E)[ 1 - / i (S)][e- < ( A '- A S +B(A "- A S)) - 1]} ^ 
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where 


T(E) = Tt{T l (E) G r (E) T r {E) G a (E)} ( 6 ) 

is the Landauer transmission coefficient at energy E. 
Here G r ( a \E) are the retarded (advanced) projections 
of the system Green function in absence of the counting 
fields, and Tr(E) is the electron dissipation matrix at 
energy E due to coupling to contact K (K = L,R). The 
size of the matrix is that of the molecular subspace of 
the problem. Below we consider systems for which ex¬ 
pression (0 is satisfied. In these systems, the electron 
dissipation rate matrices are always diagonal in the local 
basis. We denote by T^ and T# the parameters charac¬ 
terizing the electron escape rates into the left and right 
contact, respectively. 

The particle and energy average currents and noises 
can be directly obtained from the adiabatic potential U, 
Eq. 0, as 



where K = L,R and M = P, E. Explicitly, the average 
currents read 


I M = 7f = -I" 


(9) 


dE 
2i r 


OmT(E) [f L (E) - f R (E)] 


while the noises read 


qMi M2 0M1 M2 qMi M2 0M1M2 qM± M2 

° = D LL — D RR — ~ D LR — —°RL 


( 10 ) 


_ qM 1 M 2 1 nMiM 2 
shot ' ^therm ’ 


W = (hl — hr)I p , and the average heat taken from the 
hot reservoir which fuels the device, Q = — (I E — / 1 rI p ), 
namely fj = W/Q. It is upper bounded by the Carnot 
efficiency fj < 1 — Tl/Tr. The fluctuating efficiency on 
the other hand is defined as the ratio between the fluctu¬ 
ation power w/t and heat flow q/t measured at the level 
of a single experiential realization of duration t, namely 
rj = w/q. Efficiency fluctuations are not bounded and 
are characterized by the rate J (rj) at which the probabil¬ 
ity to observe a given efficiency 77 decays during a long 
measurement realization [48|, |491 | 

^) *=°°exp{-J(r))f}. (13) 

This rate is called the large deviation function (LDF) of 
efficiency. It can be derived from the heat and work FCS 
obtained from the energy and heat FCS 0 as follows. 
The heat entering the system from the hot (cold) reser¬ 
voir is given by the right (left) energy current minus hr 
(hl) times the right (left) particle current. At steady- 
state, the particle and energy currents are equal (but 
with opposite signs) at the two interfaces. Therefore, by 
the first law of thermodynamics, the work generated by 
the particles moving across the system is equal to the 
sum of the heat from the left and right reservoir which 
is thus hr — hl multiplied by the right particle current. 
This means that if A q counts the heat from the hot reser¬ 
voir and if A w counts the work, we get that the heat and 
work FCS reads 

U = iJ ^ln(l + T(E) (14) 

x |/ l (F)[ 1 - f R (E)}[e- i{[E - ,1R]XQ - [wwd A w) _ !] 
+f R (E)[ 1 - f L (Ej\[e +i{[E -^ R]XQ -^ L -^ R]Xw) - 1]} J . 


where the shot and the thermal (equilibrium) noise re¬ 
spectively read 

^f 2 =/ ^0 „iOm 2 (11) 

X T(E) (1 -T(E))(f L (E)-f R (E)] 2 

e= e / d Z° Ml o M , ( 12 ) 

K=L,R J 

x T(E) f K (E) [1 -f K (E)]. 

Expressions m-m are exact for any non-interacting 
system bi-linearly coupled to two contacts. 

III. EFFICIENCY FLUCTUATION 

In order to operate as a thermoelectric junction the 
small quantum system is embedded between two leads L 
and R with T L < Tr and hl > dR ( see Fig- CO 1 )- The 
macroscopic efficiency of such a junction is defined as the 
ratio between the average power generated by the device, 


Introducing the slightly modified version of the adiabatic 
potential, <j> = —iU , and redefining the counting fields as 
7 = iXw and A = iAq, we get that 

0(7, A ) = / ln ( 1 + T 

x {f L (E)[ 1 - f R (E)][e-U E ~^ ]A-[^-*«b) _ i] ( 15 ) 
+f R (E)[ 1 - f L (E)][e + ^ E -^ x -^ L ~^ - 1]} 

Note that the fluctuation theorem symmetry holds 

■a(7,a)=4.(-T- 7 ,_L-_L-a) (16) 

as can be verified using the property 

f R (E)[ 1 - / L (F)]e( [B_MRl( ^ _ ^” A)_[MI ' _wll(_, t _7) ) 

= f R (E)[ 1 - f L (E)\e — + ^TT e -([E-Mfl]A-[ M L-MBb) 

(17) 

= f L (E)[ 1 - f R (E)\e ~ ([E - flR]x - [,iL - ,1Rh) . 
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The efficiency LDF is finally obtain by setting A = rj_ 
and minimizing <j> relative to the field 7 , namely |48l |4S 


J(rj) = — min 0 ( 7 , 777 ). 
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(18) 


The convexity of (fl5l) together with the fluctuation theo¬ 
rem symmetry has been used in classical systems to 
prove two important results. First, the single minimum 
in J(rj) (i.e. the most probable efficiency) corresponds 
to the macroscopic efficiency 7 , second, the single maxi¬ 
mum in J(rf) (i.e. the least likely efficiency) corresponds 
to the Carnot efficiency 1 — T L /T R [H, |49|. By showing 
that the fluctuation theorem symmetry (fl6l) holds for the 
adiabatic potential of quantum junctions, we thus gener¬ 
alized these remarkable results to the quantum realm. 

In the limit of weak system-lead coupling, T = + 

r# —> 0, Eq. m reduces to the QME approach predic¬ 
tion [7ll l9~3j 


0(7, A) = 


(19) 


e(- £ # 1 + 


m 2 


r^E s )rf(E s ) 


x [f L {E s )[ 1 - f R (E s )][e-U E °-^ x -^-»*™ - 1] 
+f R (E a )[ 1 - f L (E s )][e + U E °-^ x -^™™ - 1]} 



Here ... is the sum over the eigenorbitals of the 
system with eigenenergies E s , and r s (E s ) = Tf (E s ) + 
r *(e s ) is the total escape rate from the eigenorbital s 
evaluated at energy of the the orbital. The quasi-classical 
nature of this result is manifest since Eq. (IT51) disregards 
the reservoir induced correlations between the eigenor¬ 
bitals of the system. This form of adiabatic potential 
was used in Refs. [H, [49| together with (TTH1) to calculate 
efficiency fluctuations in a photoelectric device. 


IV. NUMERICAL EXAMPLES 

We now compare the efficiency fluctuations (fT51) pre¬ 
dicted using the NEGF heat and work FCS m with 
the QME prediction (full) . Since we exclusively consider 
non-interacting models, we emphasize that the NEGF 
treatment is exact while the QME approach is an approx¬ 
imate approach only valid in the weak coupling limit to 
the contact and which neglects coherences between sys¬ 
tem eigenstates (we use the rotating wave approxima¬ 
tion to guarantee positivity). The discrepancies between 
these two approaches can thus be attributed to broaden¬ 
ing effects induced by strong coupling and to eigenbasis 
coherences. 

The calculations are performed by numerically evalu¬ 
ating the adiabatic potential <(>( 7 , 777 ) (using Eo. flTSl) for 
the NEGF and Ea. (fl9l) for the QME) and numerically 
minimizing it as a function of the counting field 7 for a 
fixed value of the efficiency 77 according to Eo. lfTSl) . 
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FIG. 2: (Color online) Efficiency LDF for a two-level bridge 
calculated within the NEGF (solid line, blue) and the QME 
(dashed line, red) approaches. See text for parameters. 


Unless specified otherwise, the parameters of the cal¬ 
culations are Tl = 100 K, T R = 600 K , /xl = 0.02 eV 
and /ijx = 0. We use the wide band approximation which 
assumes that the electron escape rates and are 
energy-independent constants. The NEGF calculations 
were performed on an energy grid spanning the region 
from — 1 to 1 eV with step 10 -5 eV. 

We start by considering the two-level bridge model 
depicted in inset in Fig. [2] when the system is weakly 
coupled to the contacts. The position of the levels is 
e 1 = £2 = 0.1 eV, the electron hopping parameter is 
t = 0.05 eV, and the electron escape rates are = T# = 
2 • 10 -4 eV. As expected, in this regime both the NEGF 
and the QME predictions for the efficiency fluctuation co¬ 
incide (compare the solid and the dashed lines in Fig. [2]). 
Large values of J{rf) indicates unlikely efficiency fluctu¬ 
ations while the minimum is the most likely efficiency 
77 corresponding to the macroscopic efficiency considered 
in traditional thermodynamics. Although hardly seen on 
this figure, the most unlikely efficiency is located at the 
Carnot efficiency 1 — Tl/T r ss 0.83. The probability dis¬ 
tribution in this regime is thus quite narrowly centered 
around the most likely efficiency. 

We consider two types of junctions, a two-level bridge 
(top inset in Fig. [3]) and a single level junction coupled to 
an isolated orbital (bottom inset in Fig. [3]). Both junc¬ 
tions are in regimes where the system is strongly coupled 
to the contacts. The latter is the simplest model often 
used to describe destructive interference effect in trans¬ 
port through a junction (see e.g. Ref. Izi). The position 
of the levels is £\ = £2 = 0.12 eV, the electron hopping 
parameter is t = 0.05 eV, and the electron escape rates 
are Tl = V R = 0.1 eV. Figure [3ji shows that QME re¬ 
sults of the two models are identical. This result stems 
from the fact that in the rotating wave approximation, 
the QME neglects coherences in the system eigenbasis 
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FIG. 3: (Color online) Efficiency LDF for a two-level bridge 
(top inset; solid line, blue) and a single level coupled to an 
isolated orbital (bottom inset; dashed line, red), calculated 
within the (a) QME and (b) NEGF approaches. The vertical 
dashed line shows the Carnot efficiency. See text for param¬ 
eters. 


[7bl ;86]. Fig. [3 Jd shows the exact efficiency fluctuations 
for the two models. The interference effects responsible 
for the discrepancy between the two curves do not signif¬ 
icantly alter the qualitative shape of the efficiency LDF. 
However, when comparing Figs. 0t and b, we note that 
the broadening effects resulting from the strong coupling 
to the contacts clearly tend to increase the magnitude of 
the efficiency fluctuations and also intensifies the asym¬ 
metry of the fluctuations around the most likely value. 
We note that even the most likely value is affected. The 
least likely value is nevertheless still exactly located at 
the Carnot efficiency. 

We now turn to the donor-bridge-acceptor (DBA) 
junction depicted in the inset of Fig. [I] This setup en¬ 
ables to study the effect of intra-molecular interference 
on efficiency fluctuations. We see that the trend pre¬ 
dicted by the QME, when moving from constructive to 
destructive interference (solid to dashed to dotted line), 
is the opposite of the real trend obtained using the exact 
NEGF. It is interesting to observe that destructive inter¬ 
ference tend to increase the most likely efficiency but at 
the same time significantly increase the magnitude of the 
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FIG. 4: (Color online) Efficiency LDF for a donor-bridge- 
acceptor junction calculated within the (a) QME and (b) 
NEGF approaches. Results are shown for constructive in¬ 
terference (s = t; solid line, blue), single path (s = 0; dashed 
line, green), and destructive interference (s = — 0.81; dotted 
line, red). The vertical dashed line shows the Carnot effi¬ 
ciency. Other parameters are as in Fig. [3] 


efficiency fluctuations. In other words, the performance 
of the junction increases but at the cost of becoming less 
reproducible. 

As a final example we consider a single level junc¬ 
tion (see inset in Fig. 0. Within the QME approach, 
the efficiency does not fluctuate in this model because 
heat and work are directly proportional to each other, 
a condition known as tight coupling [50]. However, the 
NEGF approach breaks the tight coupling condition due 
to the hybridization of the molecular level with the states 
in the contacts. The position of the level is taken as 
e = 0.1 eV and Figure [5] shows the results of calcula¬ 
tions for several strengths of the system-reservoir cou¬ 
pling: V L = r R = 0.1 eV (solid line), 0.05 eV (dashed 
line), 0.01 eV (dash-dotted line), and 0.001 eV (dotted 
line). As T —> 0 (weak coupling limit) the distribution be¬ 
comes very narrow and centered around the macroscopic 
efficiency (/x L - fi R )/ (e - fi R ). 

We now discuss ways to relate the efficiency LDF to ex¬ 
perimentally measurable characteristics of the junction. 
For the setup sketched in Fig. [T|r, the average power and 
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FIG. 5: (Color online) Efficiency LDF calculated within the 
NEGF for a single level junction. The results are shown 
for several level-contacts coupling strengths ranging from the 
strongest (solid line, blue) to the weakest (dotted line, black). 
The vertical dashed line shows the Carnot efficiency. See text 
for parameters. 

the heat flux from the hot reservoir are 

W=AfiI p (20) 

Q = - {l E - HR I P ) , (21) 

where I p and I E are defined in Eq. © and A/x = hr — 
[ir. In the linear response regime (obtained by lineariz¬ 
ing the Fermi distributions in 1/Tr^ and Hl(r)/Tr(r) 


around equilibrium hl — Hr = Ep and Tr = Tr 

= T), 

we get that 




W G An 2 + L Afi A/3 

(22) 


Q « R A/jl + F A/3 

(23) 

where A/3 = 

1/Tl - 1 /Tr and 


G 

-/gmmT 

(24) 

L 

/* d L' 

= J ^ T(E)f(E ) (E — hr) 

(25) 

R 

= ] 27i T( E )f'(E) T »* 

(26) 

F 

= - j T(E) f'(E) (E - hr) 2 - 

(27) 

Here f'(E) 

= [d/dx l/(e x + 1)\ x =(e-e f )/t and 

R = 


L/Tr. The coefficients in (1!H1) are related to experimen¬ 
tally measurable quantities. Indeed, G is the electrical 
conductance, and if k denotes the heat conductance and 
S the Seebeck coefficient, we have that 

T l Tr ; S= GT l Tr (28) 



FIG. 6 : (Color online) Efficiency LDF for a single level junc¬ 
tion (see inset in Fig.0 calculated for experimentally relevant 
parameters. The predictions of the exact NEGF calculations 
(solid line, blue) are compared to the linear response pre¬ 
dictions (1291 (dashed line, black), and to the Gaussian ap¬ 
proximation predictions ( 1331 ) (dotted line, red). See text for 
parameters. 

Thus following Ref. [If!, in the linear response regime 
the efficiency LDF can be expressed in terms of these 
measurable quantities as 

_ [t](k AT + G STr Afi) + G S AT A/x + G Ap 2 ] 2 

J(??) ~ 4 [r, 2 rT l Tr + 2 V GST l T r A/i + GTl A p 2 ] ’ 

(29) 

where AT = Tr — Tr . 

We now attempt to estimate the efficiency LDF be¬ 
yond the linear regime, solely in terms of the particle and 
energy nonequilibrium currents and noises, Eqs. ®-C2j). 
Note that in molecular junctions, the particle and energy 
currents as well as the particle noise are experimentally 
measurable [871489 L 91 and the energy noise will soon be¬ 
come measurable 94 isi- To do so, we approximate the 
cumulant generating function (fl5l) by a quadratic expan¬ 
sion in counting fields 7 and A around point 7 = A = 0. 
This is a Gaussian assumption which leads to 

0(7,777) ~ «7 2 + &7, (30) 

where we used the fact that </>( 0 , 0 ) = 0 , and defined 

a S** + ^-l‘f-’Af s PP (31) 

- V {Hl ~ Hr[^ ~ v\) S PE 
b = -r]I E + (ii l - hr[1 - rj])I p , (32) 

which are solely expressed in terms of the measurable 
nonequilibrium particle and energy fluxes, Eq. ®, and 
of the nonequilibrium noise characteristics of the junction 
m■ Within this Gaussian approximation, we find that 


k = 


(33) 
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We have thus shown that efficiency fluctuations are 
experimentally measurable close to equilibrium (1291) and 
in the Gaussian approximation <1331) . We now verify 
the validity of these approximations in Fig. [G] where the 
exact efficiency LDF, Eqs. m and m, is compared 
to the linear response, Eq. (l29l) . and the Gaussian ap¬ 
proximation, Eq. ®, result for an experimentally rel¬ 
evant set of parameters: Tl = 300 K, Tr = 350 K, 
Hl = 0.002 eV, /is = 0 We also set e = 0.1 eV and 
r L = = 0.1 eV. We see that near the minimum corre¬ 

sponding to the macroscopic efficiency, the three curves 
coincide (thus justifying the use of linear response to es¬ 
timates of average quantities). At the same time, the ef¬ 
ficiency fluctuations are poorly captured by the linear re¬ 
sponse approximation (dashed line) but reproduced quite 
well by the Gaussian approximation (dotted line). 

V. CONCLUSION 

We studied the thermoelectric properties of nanoscale 
junctions. Since stochasticity and quantum coherence are 
expected to be important at small scale, we proposed to 
characterize the performance of these devices by study¬ 
ing efficiency fluctuations rather than the widely used 
figure of merit which is intrinsically restricted to the lin¬ 
ear regime. We provided a systematic procedure to study 
efficiency fluctuations which accounts for all quantum ef¬ 
fects and is based on the work and heat FCS obtained 
within the NEGF formalism. As predicted for classical 


dynamics in Ref. [48|, the most likely efficiency coincides 
with the macroscopic efficiency, while the least likely effi¬ 
ciency corresponds to the Carnot efficiency. We used sim¬ 
ple models with realistic molecular junction parameters 
to compare our NEGF based method to the commonly 
used QME approach. We showed that the latter may 
fail qualitatively for strong system-reservoir coupling due 
to its inability to properly account for quantum coher¬ 
ences in the system. We finally proposed a method to 
estimate efficiency fluctuations using the experimentally 
measurable particle and energy nonequilibrium currents 
and noises. Linear response and Gaussian approxima¬ 
tions were proposed as ways to construct efficiency fluc¬ 
tuations from experimental measurements. We showed 
that while linear response approach, often used in the 
experimental literature to discuss thermoelectric proper¬ 
ties of junctions, captures the macroscopic efficiency, it 
fails to account for the efficiency fluctuations. At the 
same time, the Gaussian approximation was shown to 
work very well within experimentally relevant range of 
parameters. 
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Appendix A: Cumulant generating function of a multilevel non-interacting system 


Here we derive the general form of the adiabatic potential U (A) for a noninteracting n -level system. For simplicity, 
we consider the specific case of one particle counting field A in the left molecule-contact interface. Multiple counting 
fields and/or energy FCS are formulated similarly. We first find expression for the potential derivative in the counting 
field, Eqs. (0 and flxj), and then integrate it in the field to get the potential itself. 

We start by writing the dressed Green Function G(A) as a 2n x 2 n dimensional block matrix in the Keldysh contour 


roc c< 

C(A, =U 4 


(Al) 


inverse of which is 


G~\ A) = 


-iV L (f L (E) - 1/2) - ir R (f R (E) - 1/2) + IE-H m ie iX V L f L (E) + iV R f R {E) 

-ie~ iX r L (l - f L (E)) - iT R {l - f R (E)) -iT L (f L (E) - 1/2) - iV R {f R {E) - 1/2) — IE + H M 


We will use Jacobi’s formula for the derivative of the determinant of a matrix that in our case reads 


d_ 

d\ 


det(G" 1 (A)) = Tr |adj (G^A)) ^G" X (A) j 


(A2) 

(A3) 


where adj(M) denotes the adjugate matrix of a matrix M (adj(M)M = I det(M) = M adj(M)). For our consideration, 
it will be important to work with special submatrices of G _1 . For an n x n matrix M, we define to be the 

(n — 1) x (n — 1) matrix that is obtained from M by removing the jth row and the ith column. In this notation 







the (i,j)-matrix element for the adjugate of M can be expressed as adj(M)^ = (—1) J+J det(M(j|i)). Also below 
M\j i .. .j r \i i... i r ] will denote the submatrix of M composed of rows j i... j r and columns i\... i r . 

The first step in the derivation is to obtain from Eq. (1A2I) 


dX 


G~ 1 {X) 


and utilizing Eq. (IA3I) calculate 


0 -e iA T L / L (U) 

~e~ ix T L {l -f L {E)) 0 


(A4) 




Tr J 

adj (G 1 (A)) 


0 -e iX Y L f L {E) 

1 

l 

det(G" 1 (A)) 

—e 

~ lX r L {l - f L (E)) 0 

J 

Tr / 

r g%(e) g<(e) i 

0 ^e lX T L f L {E)] 

l 

[g>(£) Gi(E)\ 

[-e- lA T L (l -f L (E)) 0 


=Tr {G<(U)(-e- jA )r L (l - f L (E)) + G>(E)(-e iX )T L f L (E)} 


=H X l{E). 


(A5) 


Using this last result in Eq.Q and integrating with respect to the counting field A leads to 


f dE ln 

det(G _1 (A))" 

/ 2t r 

det(G _1 (0))_ 


(A6) 


where we used the known property £7(0) =0. Eq. (IA6I) is the first important result. 

We now have to evaluate the determinants inside the logarithm in Eq. (IA6I) . The determinants can be evaluated 
after applying elementary transformations to G _1 . First, we notice that we can write 


det(G _1 (A)) 


-E>(£) + G a,_1 -*( 1 - e iX )Y L f L {E) + E <(E) 

*( 1 - e~ lX )T L [l - f L (E)) + E >(E) -G r ~ 1 - E >(E) 


(A7) 


where G’’ 1-1 = IE — Hm + i(T l + Y R )/2 and G a,_1 = (G r,_1 )t, by adding and subtracting to each submatrix in Eq. 
(IA2D appropriate matrices. Then we add to the *-th row, i < n the (n + i)th row of the matrix. After which, on the 
resulting matrix, we add the (n + j)-th column to the j -th column for each j < n. This leads to 


det(G X (A)) 


i( 1 - e~ iX )V L {l - f L (E)) + G a,_1 -i( 1 - e iX )T L f L (E) + i( 1 - e~ iX )(l - f L {E))T L 
i( 1 - e~ iX )r L {l - f L (E)) + E >(E) i( 1 - e~ iX )r L (l - f L (E )) - G r ~ 1 


Setting A = 0 we arrive at the result for the first of the determinants in Eq. (IA6I) 


det(G _1 (0)) = det(G a,_1 ) det(—G r ’ _1 ), 


(AS) 


(A9) 


To get the second determinant in Eq. (IA6I) we have to work with the general form of Eq. (IA71) . Explicit evaluations 
lead to an expression that can be grouped in powers of (1 — e ~ lX ) and (1 — e lX ) of at most n power. In particular, 
noticing that (1 — e _lA )(l — e lX ) = (1 — e~ lX ) + (1 — e lX ) we can write 


det(G _1 (A)) = det(G a,_1 ) det(-G r ’- 1 ) + ^(1 - e" a ) s a s + (1 - e iX ) s b s , 

s=1 


where a s and b s are the coefficients of the polynomial given by 


97 


a s =(ini-f L (E)Y J2 E (-ir (Q+n)+ffW) det(r i [a|/3])det(7V(a + n|/3)) 

aeQs,n P£Qs,n 

bs =(-i) s {h(E)) s J2 E (—1)' t( “ )+<t(/3+ ” ) det(Tz,[a|/3]) det(M(a|/3 + n)) 

a&Qs,n P£Q s ,n 

where N and M are 2 n x 2 n matrices given by 


E>(E) + G a,_1 

iT R f R (E) 

, M = 

'-E >{E)+G a ’~ 1 E<(E) 

E>(E) 

-G - E >(E) 

-*r«(l -/«(£?)) —G r, ~ 1 — E>(U)_ 


(A10) 


(All) 

(A12) 


(A13) 
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Q s , n is the set of s-tuples (i 1 , ■ ■ ■, i s ) of natural numbers such that 1 < i\ < 12 < • ■ ■ < i s - 1 < i s < n, <j(a) = on 
for a € Q s ,n > an d a + n = (aq + n,a 2 + n,... ,a s + ri). From equations (1A11I) and (1A12I) we have in particular 
a n = ((1 ^ f L (E))f R (E)r det(T L T R ) and b n = ((1 - f R (E))f L (E)) n det(T L T R ). Eqs. (|A10j) - (|A12| ) give the most 
general form for the second determinant in Eq. (IA6I) . 

Finally, substituting Eqs. (IA9I) and (IA10I) into (IA6I) we obtain the general form for the adiabatic potential 


U{\) = i 


dE 
27r 


In 




(1 -e~ iX y 


—J det(G a ’ _1 ) det(— G r - *) 


b s (1 - e iX ) s 
det(G a > _1 ) det(—G r ’ _1 ) 


(A14) 


with the coefficients a s and b s given by Eqs. (IA11I) and (IA12I) . 

We now consider two specific examples where we can recover the expression derived in Ref. l69t for a single level 
junction from our general result, Eq. (IA14I) . 


1. One level. 

By direct computation of Eqs. (IA11I) and (1A12I) we find 


O! = r L r R (i - f L (E))f R (E) 61 = r L r*(i - me))Me). 


Thus 


det(G _1 (A)) 

det(G _1 (0)) 


=1 + Ga.-^gr.-i) ((! - e'^X 1 - h(E))ME) + (1 - e iA )(l - f R (E))f L {E)) 

=1 + G r (E)T R G a (E)T L ((e a - 1)(1 - f L (E))f R (E) + (e~ lX - 1)(1 - f L (E))f R (E)] 


which yields the expression for adiabatic potential derived in Ref. 


2. n-level system coupled to the contacts through single orbitals. 

In this case and are nxn matrices with all entries equal zero but one element in the diagonal. Examples 
of systems of this kind are the two level bridge (see inset in Fig. EJ> or D-B-A type of the junction (see inset in 
Fig. [4]). Here we can take [r i]ij = SijSn'yL and [Tr]^- = S n j6i n "f R , which results in a s = b s = 0 for s > 1 and 


ai =i( 1 - f L {E)){- l) 1+n+1 det(r L [l|l]) det(AT(l + n|l)) 


=i(l-f L (E))(-ir lL 


-S>(|1)( J E7) -F G^’-Hll) 
-*Fr(1|1)(1 - f R (E)) 


iT R f R (E) 

—G r ’“ 1 (l|) — E>(1|)(I?) 


=*(i - f L (E))(-ir^ L (-ir +2n -H 1R f R (E) 


—E > (n|l)(£ : ) + G a,-1 (n|l) 0 

-iF«(l|l)(l - f R (E)) —G r,_1 (l|n) - E>(1| n)(E) 


=(1 - fL{E))f R (E)'yL n /R det(G a ’ _1 (n|l)) det(—G r ’ _1 (l|?r)) 


(A15) 


where we used E- 5 " (rz| 1) (E) = E > (l|n)(Fi) = 0. Similarly 

h =h(E)(l -/ fl (F;))7 L 7Rdet(G a ’- 1 (l|n))det(-G"’- 1 (77|l)). (A16) 

From the definition of the adjugate we have (—1) 1+ ” det(G a,_1 (l|n)) = adj(G a,_1 ) n i, 
(—1)” +1 det(—G r, ~ 1 (n|l)) = adj(—G r,_1 )i n . Also in this particular case adj(G° ,_1 )i„ = adj(G r, ~ 1 )i„ 
and adj(G a ’~ 1 )„i = adj(G r, ~ 1 ) ri i. Finally, substituting Eqs. (I A 1511 and (IA16I) into Eq. HAITI) and rearranging 
terms, we get 


U{ A) =i 



l + Ti-{G r (E)T R G a (E)T L } 


x ((e iA - 1)(1 - f L (E))ME ) + (e~ iX - 1)(1 - f L {E))f R (E )) 


(A17) 


Eq. (IA17I) is the Levitov-Lesovik formula for a multi-level system. 
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